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Abstract 

It is argued that to arrive at a quantitative description of the surface tension of a 



liquid drop as a function of its inverse radius, it is necessary to include the bending 
rigidity k and Gaussian rigidity k in its description. New formulas for k and k in 



the context of density functional theory with a non-local, integral expression for the 
interaction between molecules are presented. These expressions are used to investigate 
the influence of the choice of Gibbs dividing surface and it is shown that for a one- 
component system, the equimolar surface has a special status in the sense that both 
k and k are then the least sensitive to a change in the location of the dividing surface. 
Furthermore, the equimolar value for k corresponds to its maximum value and the 

in 

a short-ranged interaction potential between molecules, shows that k is negative with 

O 

a value around minus 0.5-1.0 k-^T and that k is positive with a value which is a bit 



equimolar value for k corresponds to its minimum value. An explicit evaluation using 



more than half the magnitude of k. Finally, for dispersion forces between molecules, 
we show that a term proportional to \og(R)/R 2 replaces the rigidity constants and 
we determine the (universal) proportionality constants. 



I. INTRODUCTION 

The surface tension of a simple drop of liquid has captured the imagination of sci- 
entists dating back to the pioneering work of J. Williard Gibbs [1]. This interest 
continues with the main focus of attention directed towards the description of the 
deviation of the surface tension from its planar value when the radius of the liquid 
droplet becomes smaller. Such a deviation is especially important in the theoretical 



description of nucleation phenomena [2]. The homogeneous nucleation of a liquid 
from a supersaturated vapour follows via the formation of small liquid droplets and 
the nucleation time and energy depend sensitively on the precise value of the droplet 's 
surface tension. 

A key quantity in quantifying the extent by which the surface tension of a liquid 
drop deviates from its planar value is the Tolman length introduced by Tolman in 
1949 [3]. It can be defined in two equivalent ways. In the first way, one considers the 
radial dependence of the surface tension of a (spherical) liquid droplet defined as the 
excess grand free energy per unit area: 

n = -p e V e -p v V v + a s (R)A. (1) 

When the radius R of the droplet is large, the surface tension may be expanded in 
the inverse radius: 

*.{R) = *~ + ..., (2) 

where a is the surface tension of the planar interface and where the leading order 
correction defines the Tolman length 5. In the second route to define the Tolman 
length, one considers the pressure difference Ap=pg — p v between the pressure of the 
liquid inside and the pressure of the vapour outside the droplet. For large radii of 
curvature, Ap is expanded in 1/R: 

2a 25a 

Ap -n~w + -- (3) 

The first term on the right hand side is the familiar Laplace equation [41 with the 

n 

leading order correction giving Tolman's original definition of the Tolman length [3|. 
It is important to note that this correction only takes on the form in Eq.([n]) when the 



ii 



equimolar radius |l[ is taken as the radius of the liquid drop, i.e. R = R e . Furthermore, 
with this choice of the (Gibbs) dividing surface, terms of order 0(1/R 3 ) are absent 
and the dots represent terms of order 0(1/ R 4 ). When the location of the droplet 
radius is chosen away from the equimolar radius, the Tolman length correction to the 
Laplace equation has a form different than that shown in Eq.([3]). For instance, the 
radius corresponding to the so-called surface of tension (R = R S ) is defined such that 
Eq.([3]) appears as Ap = 2a(R s )/R s . 

The determination of the value of the Tolman length for a simple drop of liquid 
has proved to be not without controversy (recent reviews are given in refs. [5|, la]). 



This is mainly due to two reasons: first, one of the first microscopic expressions for 
the Tolman length was formulated in the context of a mechanical approach which 
lead to an expression for the Tolman length in terms of the first moment of the excess 
tangential pressure profile of a planar interface [7|. However, it was pointed out by 
Henderson and Schofield in 1982 that such an expression dep ends on the form of the 
pressure tensor used and is therefore not well-defined |8l— 1 1 1 j| . Furthermore, even the 
evaluation of the Tolman length using the usual Irving-Kirkwood 12| form for the 
pressure tensor leads to incorrect results [13| and the use of the mechanical expression 
is now (mostly) abandoned. 

A second origin of controversy is simply due to the fact that for a regular liquid- 
vapour interface the Tolman length is small (a fraction of the molecular diame- 
ter), since it measures the subtle asymmetry between the liquid and vapour phase. 
Straightforward squared-gradient theory with the familiar tan/i-profile for the density 
profile, leads to a zero value of the Tolman length [14J, ll5( and it remains a challenge 
to distinguish its value from zero in computer simulations |16|-|20|. Nowadays, those 



computer simulations that have succeeded in obtaining a value different from zero 
indicate that its value is negative with its magnitude around one tenth of a molecular 



diameter 
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26| and error bars usually somewhat less than half that number. 



The sign and magnitude of the Tolman length for a regular liquid- vapour interface 
are corroborated by a large number of different versions of density functional theory 
(DFT), which has proved to be an invaluable tool in the theoretical description of 



inhomogeneous systems [27|-|30j. Quite surprisingly, the details of the density func- 
tional theory at hand do not seem to matter that much p, 13JJ] and one ubiquitously 
finds that the Tolman length is negative with a magnitude comparable to that ob- 
tained in simulations. This includes results for the Tolman length from van der Waals 
squared-gradient theory [32|, |33|, density functional theory with a non-local, integral 



expression for the interaction between molecules (D 
tional theory with weighted densities (DFT-WDA) 



T-LDA) 



lalaikz 



371 ] , density func- 



using Rosenfeld's 
(DFT-FMT) [23 



m 



3jJ and density functional theory 
38 1 fundamental measure theory for the hard-sphere free energy 



39]. 



All in all, there now seems to be the same level of agreement between simulations 
and DFT for the Tolman length as it exists for the surface tension, with the exception 



of one particular type of simulation result. In refs. 
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- |26J the Tolman length is 



determined in computer simulations of liquid droplets for various (large) radii of 
curvature, but in a different set of simulations the Tolman length is extracted from 
computer simulations of a planar interface 40|, l4l| . using a virial expression for the 



Tolman length 42j. The simulations of the planar interface lead to a Tolman length 
that has the same order of magnitude as the simulations of the liquid droplets but 
now with the opposite sign. It has been suggested that, since the interfacial area is 
much larger in the simulations of the planar interface, the presence of capillary waves 
might play an important role 21] . However, it is difficult to imagine that this would 
change the sign of the Tolman length so that the resolution to this problem remains 
uncertain. 

Another feature that ubiquitously results from the computer simulations and DFT 
calculations of liquid droplets is that the surface tension is not monotonous as a 

n 

function of the (inverse) radius (for a recent review, see ref. |6j). A maximum in the 
surface tension of a liquid droplet occurs which suggests that the surface tension is 
qualitatively better approximated by a parabola rather than by a straight line with its 
slope given by the Tolman length. This means that one needs to include higher order 
terms, going beyond the level of the Tolman length, in the expansion of the surface 
tension in Eq.(j2J). Such an expansion was first provided in the ground-breaking work 
by Helfrich in 1973 [43]. The form for the free energy suggested by Helfrich is the 
most general form for the surface free energy of an isotropic surface expanded to 
second order in the surface's curvature [43] : 

n H = fdA[a-5a.J+- J 2 + kK +...], (4) 

where J—I/R1 + I/R2 is the total curvature, K = 1/(R]_R2) is the Gaussian curvature 
and Ri, R2 are the principal radii of curvature at a certain point on the surface. 
The expansion defines four curvature coefficients: a, the surface tension of the planar 
interface, 5, the Tolman length [3J], k, the bending rigidity, and k, the rigidity constant 
associated with Gaussian curvature. The original expression proposed by Helfrich 



43] features the radius of spontaneous curvature i? as the linear curvature term 
(5 a — > 2k /Rq [a, ll3]), but in honour of Tolman we stick to the notation in Eq.fjlJ). 

For surfaces for which the curvatures J and K are constant, the Helfrich free 
energy per unit area reduces to: 

n u /A = a(J,K) = a-5aJ+-J 2 + kK + ... , (5) 



which for a spherically or cylindrically shaped surface takes the form: 

25a (2k + k) . , , 

a s (R) = a - — + K R2 ' +... (sphere) (6) 

0c(R) = ° - -£■ + ^2 + • • ' (cylinder) (7) 

These expressions indicate that the second order coefficients, which express the non- 
monotonicity of the surface tension as observed in simulations and DFT calculations 
of liquid drops, are given by the combination of the rigidity constants 2k + k and 
the bending rigidity k. Our goal in this article is to provide general formulas for the 
bending rigidities k and k using density functional theory (DFT-LDA). This work 



extends previous work by us [34], by Koga and Zeng [44j, by Barrett |45j and by 



Baidakov et al. (46] . Our formulas are subsequently applied to explicitly evaluate 
the bending rigidities and it is determined how well they can be used to describe the 
surface tension of a liquid drop (or vapour bubble) . 

The expansion of the surface tension of a liquid drop to second order in 1/R has 



not been without controversy 47H49J . Two issues have played a role here. The 



first issue concerns the fact that when the interaction between molecules is suffi- 



ciently 
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ong-ranged, the expansion in 1/R may not be analytic beyond some term 



511 ]. In particular, for dispersion forces the second order correction has the 
form log(R)/R 2 rather than 1/R 2 and one could argue that the rigidity constants are 
"infinite" . Nowadays, this point is well-appreciated and no longer source of contro- 
versy. In this article we come back to this issue and provide explicit expressions for 
the second order correction to replace the expansion in Eq.([6]) or ((7|) for dispersion 
forces. 

A second issue argues that even for short-ranged interactions, which are mostly 
considered in simulations and DFT calculations, the second order term might pick 



up a logarithmic correction of the form log(R)/R 2 47H49J . The reasoning behind 
this focuses on the fact that for a spherical droplet, the second order contribution 
to the free energy, i.e. the expression in Eq.(jSJ) multiplied by the area A = Air R 2 
is independent of R, which might be an indication that it should be replaced by 
a logarithmic term. The most compelling argument against this reasoning lies in 
the fact that the same argument applied to a cylindrical interface would lead to 
the conclusion that already the linear term in 1/R (Tolman length) would pick up 
logarithmic corrections. Although the issue is not completely settled, the presence 



of a logarithmic correction for short-ranged interaction has not been observed in 
simulations or demonstrated in calculations either in mean-field theory (DFT) or in 



Statistical Mechanics J42|. Also in this article, we inspect (numerically) the possible 
presence of a logarithmic correction to the second order term in the expansion of the 
free energy of a liquid drop and find no evidence for its presence. 

Our article is organized as follows: in the next section we discuss the density 
functional theory that is considered (DFT-LDA) and use it to determine the surface 
tension a s (R) of a liquid drop and vapour bubble. In Section [TTT| the free energy is 
expanded to second order in 1/R for a spherical and cylindrical interface which al 



ows 



the formulation of new, closed expressions for the rigidity constants k and k [34], |45 |. 
An important feature addressed is the consequence of the choice made for the location 
of the dividing surface (the value of R) on the value of the bending rigidities. The 
formulas for k and k are explicitly evaluated using a cut-off and shifted Lennard- Jones 
potential for the attractive part of the interaction potential. Since the evaluation of 
these expressions requires numerical determination of the density profile, we supply in 
Section IIVI an accurate approximation based on squared-gradient theory to evaluate 
5, k and k from the parameters of the phase diagram only. In Section |V] we consider 
the full Lennard- Jones interaction potential and determine its consequences for the 
expansion of the free energy in 1/R. We end with a discussion of results. 

II. DENSITY FUNCTIONAL THEORY 

The expression for the (grand) free energy in density functional theory is based on 
the division into a hard-sphere reference system plus attractive forces describe d by 



an interaction potential £/ att (r). It is the following functional of the density r 27H30|: 



m = Jdr [ fUp) - Mr) ] + ~ Jdnjdfu U M {r) p{r x )p{r 2 ) , (8) 

where p is the chemical potential. For the free energy of the hard-sphere reference 



system /hs(p)> we take the well-known Carnahan-Starling form 52| : 



fhs (p) = k B T p ln(p) + k B T p , (9) 

where 7] = (7r/6) p d 3 with d the molecular diameter. The Euler-Lagrange equation 
that minimizes the free energy in Eq.flH]) is given by: 

V- = fL(p) + [dri2 U m (r) p(r 2 ) . (10) 



For a uniform system, the Euler-Lagrange equation becomes: 

H = fM-2a Pi (11) 

with the van der Waals parameter a explicitly expressed in terms of the interaction 
potential as 

a = - l -Jdf l2 U^{r). (12) 



Using the expression for the chemical potential in Eg. lull) , the bulk pressure is ob- 
tained from Q = —pV leading to the following eguation of state: 

fc B Tp (1+T7 + T7 2 -T? 3 ) 2 

p = (TT^ ap ■ ( 13 ) 

Next, we consider the implementation of DFT in planar and spherical geometry. 

Planar interface 

When the chemical potential is chosen such that a liguid and vapour phase coexist, 
/i = /i coex , a planar interface forms between the two phases. The density profile is then 
a function of the coordinate normal to the interface, p(f) = po(z). In planar geometry, 
the Euler-Lagrange eguation in Eg. ffTQl) becomes: 

Aicoex = /hs(Po) + df 12 U att (r) p {z 2 ) . (14) 



The surface tension of the planar interface is the surface free energy per unit area 

(a = (n + pV)/A[A}): 

oo 

a = -- Jd Zl Jdf 12 U att (r) r 2 (l - s 2 ) p' { Zl ) p' {z 2 ) , (15) 

— oo 

where z 2 — Z\ + sr and s = cos9i 2 . 

A Spherical Drop of Liquid 

When the chemical potential \x is varied to a value off- coexistence, spherically shaped 
liguid droplets in metastable eguilibrium with a bulk vapour phase may form. Such 
droplets are termed critical droplets. The radius of the liguid droplet is taken to be 
egual to the equimolar radius, R = R e (l|, which depends on the value of the chemical 
potential chosen, and is defined as: 

/Att 
dr r 2 [p s (r) - p v ] = — R 3 e ( Pi - p v ) . (16) 



The (grand) free energy for the formation of the critical droplet is given by: 

AQ _ Q+ Pv V ApR 

— = —4— = —a" + as{R) ' (17) 

with p v the vapour pressure outside the droplet and p? = p v + Ap is the liquid pressure 

inside (see the remark below, however). The surface tension of the critical droplet is 

the quantity that we wish to study and this equation provides a way to determine it 

from AQ. 

In spherical geometry, the free energy density functional in Eq.Q is given by: 



A 





jdn \jj;J [ fhs(ps) - W,(ri) ] 


00 2 

+ 2J dTl \Rj J™ f/att ( r )^( ri )^( r2 ) 



n 
with the Euler-Lagrange equation that minimizes the above free energy equal to: 



P = fL(Ps) + Jdr 12 f/ att (r) p s (r 2 ) . (19) 

The procedure to determine a s (R) as a function of R is as follows: 

(1) First, the bulk densities po,£ an d po,v an d the chemical potential at two-phase 
coexistence, p coex , are determined by solving the following set of equations: 

/ (P0,v) = /^coex , / {P0,e) — /^coex , f{P0,v) ~ Pcoex p0,v = f[.P0l) ~ Pcoex Po,£ , (20) 

where we have defined f(p) = f hs (p) — ap 2 . The bulk density difference is denoted as 
^P = Po,£-Po,v and the pressure at coexistence is simply p CO ex = -f{Po,t/v) +P>coexPo,t/v ■ 

(2) Next, the chemical potential /i is varied to a value off- coexistence. For /i>/x coex 
liquid droplets are formed (R > 0) and when p < p coex we obtain bubbles of vapour 
(R < 0). For given temperature and chemical potential p the liquid and vapour 
densities pt and p v are then determined from solving the following two equations 

f(p v )=p, f\pe)=P, (21) 

with the corresponding bulk pressures calculated from 

Pv = ~f(pv)+ ppv, Pe = -f(pe)+ PPl- (22) 

It should be remarked that far outside the droplet (r— >oo), the density (or pressure) 
is equal to that of the bulk, p s (oo) =p v , but that only for large droplets is the density 
inside the droplet (p s (r — 0)) equal to its bulk value {pi). 

8 
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FIG. 1: Phase diagram as a function of reduced temperature and density. The solid lines 
are the liquid-vapour densities at two values of the reduced LJ cut-off radius (solid circles 
indicate the location of the critical points). Square symbols are simulation results from 



ref. [53|]. 



(3) Finally, the Euler-Lagrange equation for p s {r) in Eq.( TT9l) is solved numerically 
with the boundary condition p s (oo) = p v . The resulting density profile p s {r) is inserted 
into Eq. ffl~6|) to determine the equimolar radius R = R e and into Eq. ffTS]) to determine 
AQ and thus a s (R). 

This procedure is carried out using a cut-off and shifted Lennard- Jones potential for 
the attractive part of the interaction potential: 

^Lj(V m in) - U L j(r c ) < r < r min 

U att (r) = \ U L ,(r)-U LJ (r c ) r min < r < r c (23) 





mm <-. T <-. T c 

r > r r 



where Uu(r) = 4e [ (d/r) 12 — (d/r) 6 ] and r min = 2e d. Figure [T] shows the resulting 
phase diagram as a function of reduced density p* = pd 3 and reduced temperature 
T* = k-^I '/e. The solid lines are the liquid-vapour densities for two values of the LJ 
cut-off radius; the square symbols are recent computer simulation results taken from 
ref. |53|. 

In Figure EJ we show the pressure difference multiplied by R/2a as a function of 
the reciprocal radius. The circular symbols are previous simulation results [21] that 
were used to determine the Tolman length from (minus) the slope at 1/R = (5~ - 



0.10 d 



2 1||). For comparison, we show the result of DFT calculations as the solid line, 



where we have taken the pressure at the center of the droplet as the liquid pressure. 



1.03 






1.02 



1.01 



1.00 




0.1 0.15 



FIG. 2: Pressure difference multiplied by R/2a as a function of the reciprocal equimolar 



211 ]. DFT calculations are 



radius d/R. Circular symbols are simulation results from ref. 
shown as the solid line (Ap = p(0) — p v ) and square symbols (Ap = pi — p v )- For the DFT 
calculations we have set the reduced temperature T* = 0.911297 and reduced LJ cut-off r c = 
2.5. The value for the reduced temperature is chosen such that the liquid-vapour density 



difference at coexistence matches the value in the computer simulations 21] 



The excellent agreement in Figure [2] is somewhat misleading since the corresponding 
values of the surface tension differ by as much as 50 %. As square symbols, the results 
of DFT calculations using pi from Eg. ( 122]) as the liquid pressure are plotted to show 
that the slight difference between p(0) and pt for small droplets has no consequences 
for the determination of S. 

In Figure El a typical example of the surface tension of a spherical liquid drop (and 
vapour bubble) is shown as a function of 1/R, with R the equimolar radius of the 
droplet. The symbols are the values for a s (R) calculated using DFT. The solid line is 
the parabolic approximation in Eq.(j6]) with values for the coefficients a, 5, and 2k + k 
calculated from formulas presented in the next Section. The behaviour of the surface 
tension is characterized by a positive first derivative at 1/R = 0, which indicates 
that the Tolman length is negative, and a negative second derivative which indicates 
that the combination 2k + k is also negative. It is concluded that the parabolic 
approximation gives a quantitatively accurate description for the surface tension for 
a large range of reciprocal radii. The determination of the full u s {R) is usually quite 
elaborate and it therefore seems sufficient to only determine the coefficients in the 
parabolic approximation to u s {R) as a function of 1/R. This is done in the next 



10 



0.22 


i i i i 


i i i 




*^^ 


^^T^. *"* "^ 




** or 
* or 


^□W *% 
^b*. ^ 


0.20 


s or 
* Jo 






* (B 


°\ 


o(R) 


* 2> 

• IB 


o\ 
o\ 
o\ 




/ i 


o\ 


0.18 


- / 6 


o\ 

°\ " 




' P 


\ 








0" 

o 


0.16 


i / i i i 


1 1 1 



-0.2 



d/R 0.2 



FIG. 3: Droplet surface tension (in units of k^T/d 2 ) as a function of the reciprocal equimolar 
radius d/R; vapour bubbles are formed for R<0 and liquid droplets for R>0. The solid 
line is the parabolic approximation to <r s (R) determined from the expansion in Section IIII1 
As a comparison, the parabolic approximation to the cylindrical surface tension cr c (R) is 
shown as the dashed line. We have set the reduced temperature T* = 1.0 and reduced LJ 
cut-off r c = 2.5. 

Section. 

III. CURVATURE EXPANSION 

In this section, we consider spherically and cylindrically shaped liquid droplets and 
expand the free energy and density profile systematically to second order in 1/R. An 
important feature of our analysis will be to not restrict ourselves to a particular choice 
of the dividing surface, but to instead leave the radius R unspecified. This will allow 
us to derive new, more general expressions and will allow for a new investigation of 
the consequences of varying the choice for the location of the dividing surface. 

To second order in 1/R, the expansion of the density profile of the spherical droplet 
reads: 



1 1 

Ps{r) = PoO) + ^Ps,i(z) + -j^Ps,2( Z ) + ■ 



(24) 



where z = r — R. The leading order correction to the density profile of the spherical 
interface is twice that of the cylindrical interface, so it is convenient to define pi{z) = 
Ps,i( z ) = 2p Cj i(z). We shall consider the expansion of the free energy of the spherical 
and cylindrical droplet separately. 
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Spherical interface 

The coefficients in the curvature expansion of the density are determined from the 
curvature expansion of the Euler-Lagrange equation in Eq. (1191) . The result is that 
the (planar) density profile po{z) is determined from Eq. (TT4|) and pi{z) follows from 
solving: 

2 

Pi = fL(Po) Pi(*i) + Jdr l2 tU(r) [pi(z 2 ) + y (1 - s 2 ) p' (z 2 ) } , (25) 

where pi = 2a/Ap p, ll5|. For the evaluation of the curvature coefficients it turns out 
to be sufficient to determine the density profiles po(z) and pi(z) only. 

The expansion for p s (r) is inserted into the expression for the free energy in Eq. (I18|) . 
Performing a systematic expansion to second order in 1/R, using the Euler-Lagrange 
equations in Eqs. (fT4"j) and (1251) . one ultimately obtains expressions for the curvature 
coefficients by comparing the free energy to the curvature expansion in Eq.Q. For 
the surface tension of the planar interface the result in Eq. (fT5"j) is recovered: 

1 °° 
a = - - jd Zl Jdf 12 £/ att (r) r 2 (l - s 2 ) p' ( Zl ) p' (z 2 ) . (26) 



For the Tolman length one obtains the following expression 34] 



8a = \ Jd Zl Jdf l2 f/ att (r)r 2 (l - s 2 )z x p' {z x )p' Q {z 2 ) - ^ jdz z p' {z) . (27) 

— oo — oo 

For the combination of the rigidity constants, 2k + k, we have: 

oo 

2k + k = - Jd Zl jdr 12 U m (r) r 2 (l - s 2 ) p' ( Zl ) Pl (z 2 ) (28) 

— oo 
oo 

- t \dz\ dfi2 f/ at t(r)r 2 (l - s 2 ) zf p' (zi) p' (z 2 ) 

— oo 
oo 

+ ^g Jdzxjdr 12 f/ att (r)r 4 (l - s 4 ) Po(2i)pd(z 2 ) 

— oo 

oo 

'cfe y^pi(z) +/ii2 2 po(z) + p Sj2 zp' (z) 

— oo 

where p Sj2 = —aApi/(Ap) 2 — 25a /Ap p, |l5| with Api = pi/ — p± tV . 

Cylindrical interface 

The analysis for the cylindrical interface is analogous to that of the spherical interface. 
Following the same procedure as for the spherical interface, the expressions for a and 

12 



5a in Eqs.( )26p and (|27[ ) are recovered and one obtains as an expression for the bending 
rigidity k: 

oo 

k = - Jdz x Jdr 12 U m {r) r 2 (l - s 2 ) p' G {z 1 )p 1 {z 2 ) (29) 

— oo 
oo 

+ l~ A /^i/rfr 12 f/ att (r)r 4 (l- S 2 ) 2 po(^i)Po(^) 

— oo 
oo 

+ Jdz [^zp[(z) + ^z 2 p' (z)+2p c , 2 zp' (z) 

— oo 

where /i C)2 = —a Api/(2 Ap) 2 p, |l5j . An expression for the rigidity constant associated 
with Gaussian curvature is then obtained by combining Eqs. fl28|) and ( 129]) : 

oo 

k = -- Jd Zl Jdf 12 U m (r) r 2 (l - s 2 ) z 2 p'oO^'oO^) (30) 

— oo 
oo 

-gg /^i/rfr 12 f/ att (r) r 4 (l - S 2 )(l - 5s 2 ) p'^z,) p' (z 2 ) 



+{p s ,2 - 4p c , 2 ) Jdz zp' (z) 



The expressions for k and k differ in two ways somewhat from previous expressions 
derived by us in ref. [341. First, they are rewritten in a more compact form with 



a printing error in ref. [34] corrected (as noted by Barrett 45] ). Second, these ex- 
pressions are derived without reference to a particular choice for the location of the 
dividing surface, i.e. for the location of the z = plane. This feature allows us to 
investigate the influence of the choice for the location of the dividing surface. As 
already known, the surface tension and Tolman length are independent of this choice 
but k and k do depend on it. 

Choice for the location of the dividing surface 

We first consider the density profile of the planar interface, obtained by solving the 
differential equation in Eq.( Tb£l) . to investigate the consequences of the choice for the 
location of the dividing surface for 5 and k. One may verify that when po(z) is a 
particular solution of the differential equation in Eq. (fT4"j) . then the shifted density 
profile 

p {z) — y p (z - zq) , (31) 
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is also a solution for arbitrary value of the integration constant Zq. However, since the 
expressions for 5 and k feature z (or z\) in the integrand, such a shift has consequences 
for the different contributions to 8 and k. To investigate this in more detail, we first 
place the dividing surface of the planar system at the equimolar surface, z = z e , which 
is defined such that the excess density is zero [l[: 

oo 

dz [poO) - Po,e @0e - z) - p ,v ©0 - z e )} = - Jdz (z - z e ) p' {z) = , (32) 

) — oo 

where 0(z) is the Heaviside function. When all distances to the surface are measured 
with respect to the equimolar plane, we need to replace z by z — z e in the expressions 
for S and k. For the Tolman length in Eq. ([2"Tj) we then find that: 

oo 

5a = - jdzx dr u £/ att (r) r 2 (l - s 2 ) (z x - z e ) p' (zi) p' (z 2 ) , (33) 

— oo 

where we have used Eq.f )32p . Now, to investigate the consequences of shifting the 
dividing surface away from the equimolar surface by a distance A, we replace z — > 
z — (z e + A) in the expression for the Tolman length in Eg. (1271) . One may easily verify 
that on account of the fact that fi\ = 2a/ Ap the Tolman length then again reduces to 
the expression in Eq. fl33l) which proofs that the Tolman length is independent of the 
choice for the location of the dividing surface. 

Replacing z — > z — z e in the expression for the rigidity constant associated with 
Gaussian curvature in Eq. (l3TJl) . we find that k simplifies to 

oo 
^equimolar = ~J \dz X df U U at t(r) T 2 (l - S 2 ) {z\ ~ Z e ) 2 p' Q (zi) p' Q (z 2 ) (34) 



-oo 
oo 



~ jd Zl jdr 12 C/ att (r) r 4 (l - s 2 )(l - 5s 2 ) p' Q {z x ) p' (z 2 ) . 

— oo 

Again, we may investigate the consequence of shifting the dividing surface by replacing 
z — > z — (z e + A) in the expression for k in Eq.f l30p . We then find that 

k — rCequimolar + 0" A . (35) 

This equation shows that k does depend on the choice for the location of the dividing 
surface. It also shows that k evaluated for the equimolar surface (A = 0), corresponds 
to the lowest possible value for k and is the least sensitive to a shift in the location 
of the dividing surface. 
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To address the influence of the dividing surface on the value of the bending rigidity 
k, we need to consider the properties of the density profile p\{z) as well. One may 
verify that when p x [z) is a particular solution of Eq. (|25|) then also 

px{z) — > pi(z) +ap' {z), (36) 

is a solution for arbitrary value of the integration constant a. Now, one may easily 
verify by inserting Eq.( |36l) into Eq. (129ft that k is independent of the value of the 
integration constant. This means that just like 5 and k we only need to consider the 
influence of the choice for the location of the dividing surface of the planar density 
profile po(z). For the equimolar surface, the expression for the bending rigidity in 
Eq. (l29l reduces to: 



k 



equimolar 



oo 

= - Jdz x Jdf 12 C/ att (r) r 2 (l - s 2 ) p' (z 1 )p 1 (z 2 ) (37) 

— oo 

1 °° 
+ - fd Zl Jdf 12 U att (r)r 4 (l- s 2 ) 2 Po(zi)p'o(z 2 ) 

z - z e ) p[(z) +2(z- z e ) 2 p' (z) 



+ f /,/, 



-oo 



Shifting the dividing surface by replacing z— >z — (z e + A) in the expression for k in 
Eq.(l29l. we then find that 

k = ^equimolar 0~ IA . {do) 

It is concluded that also the bending rigidity k does depend on the choice for the 
location of the dividing surface. The bending rigidity evaluated for the equimolar 
surface (A = 0), now corresponds to the largest possible value for k but it is again the 
least sensitive to a shift in the location of the dividing surface. 

The procedure to determine the curvature coefficients a, 5, k and k is now as follows. 
The planar profile po(z) is first determined from the differential equation in Eq. tTHl) 
with p £, p 0v , yU coex and p coex derived from solving the set of equations in Eq. (l20|) . 
From po{z), the location of the equimolar plane z = z e is determined from Eq. (l32l) and 
the curvature coefficients a, 5 and k are evaluated from the integrals in Eq.( r26|) . ( l33j) 
and (Ell), respectively. The constant pi is subsequently determined from pi = 2a/Ap 
which allows us to determine the bulk density values pi// v from pi// v =Pi/ f"(po,£/ v )- 
For given po(z) and pi, the differential equation for p±(z) in Eq. (l25|) is solved with the 
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FIG. 4: Surface tension a (in units of k^T/d 2 ) and Tolman length 5 (in units of d) as a 
function of reduced temperature. Circular symbols are the results of the full DFT calcu- 
lations in Eqs. (|26p and ()27p . The solid lines are the squared-gradient approximations of 



Section ITVT . Square symbols are simulation results for a from ref. 
(solid square) and ref. 24| (two open squares). 



and for 5 from ref. 



2ii 



boundary conditions pi(—oo) = pi t e and pi(oo)=pi tV . Finally, with pi(z) determined, 
k can be evaluated from the integral in Eq.f l37|) . 

This procedure is carried out (again) using the cut-off and shifted Lennard- Jones 
potential in Eq.( !23|) for the attractive part of the interaction potential. Figure H] 
shows the surface tension and Tolman length as a function of temperature. The 
circular symbols are the values for a and 5 calculated using DFT for two values of 
the LJ cut-off radius r c . The solid lines are the squared-gradient approximations in 
Section |IV] for r c = 2.5, 7.5, and oo. As square symbols, we show computer simulation 
results for a from ref. 53| , the single simulation result for S from ref. 21] (solid square) 
and results for 5 from simulations by the group of Binder 24j (open squares). 
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FIG. 5: Bending rigidity k, Gaussian rigidity k, and the combination 2k + k (in units of 
k&T) as a function of temperature. The rigidity constants are evaluated using the equimolar 
surface as the dividing surface. Circular symbols are the results of the full DFT calculations 
in Eqs. (|34p and (|57|) . The solid lines are the squared- gradient approximations of Section 



IIV1 Square symbols are simulation results from ref. 



24] 
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In Figured, the bending rigidity k, Gaussian rigidity k, and the combination 2k+k 
are shown as a function of temperature. The rigidity constants are evaluated using 
the equimolar surface for the location of the dividing surface. The circular symbols 
are the values for k and k calculated using DFT for two values of the reduced L J cut- 
off radius r c = 2.5 and 7.5, with the solid lines the corresponding squared-gradient 
approximations determined in the next Section. Also shown are simulations results by 
the group of Binder [24!]. Although a detailed comparison of the DFT and simulation 
results is not really appropriate due to a difference in cut-off used, the agreement in 
sign and order of magnitude is rather satisfactory. 

IV. SQUARED-GRADIENT EXPRESSIONS 

The evaluation of 5, k and k requires the full numerical evaluation of the density 
profiles po(z) and p\(z) from the differential equations in Eqs. f)14l) and (1251) . This 
procedure is quite elaborate, prompting a need for simple formulas that provide (ap- 
proximate) numbers for the various coefficients. In this section we provide a rather 
accurate approximation scheme based on the squared-gradient approximation which 
only requires the calculation of the phase diagram as input. 

The squared-gradient theory for surfaces dates back to the work of van der Waals in 



1893 54]. Its free energy functional is derived from Eq.flH]) by assuming that gradients 



in the density are small so that p(r*2) may be expanded around p(r*i). This leads to: 

Q[p] = fdf \m | Vp(f)| 2 + f(p) - pp{r)\ , (39) 



where the squared-gradient coefficient m is given by 

m = -— Jdr 12 r 2 U att (r) . (40) 

Expressions for the curvature coefficients in squared-gradient theory were formulated 
some time ago. For the surface tension of the planar interface, we have the familiar 
expression given by van der Waals [f 



a = 2m dz p' (z) 2 . (41) 

— oo 

For the Tolman length, Fisher and Wortis derived the following expression |l4J: 

oo 

5a = -2 m fdz {z - z e ) p' (z) 2 . (42) 

— oo 
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For the bending and Gaussian rigidity, one has [151 ]: 

oo oo 

k = -m J dz po{z) p[(z) + J dz -^-z p[(z) + y z 2 p' (z) + 2p c ^zp' (z) 

— oo — oo 

oo oo 

k = 2m dz z 2 p' (z) 2 + (p s>2 - 4/i Cj2 ) dz z p (z) , (43) 

— oo — oo 

which, evaluated using the equimolar surface for the location of the dividing surface, 
reduce to: 

oo oo 

^equimolar = -TTl I dz p (z) p[(z) + -j J dz (z - Z e ) p\(z) + 2 (z - Z e ) 2 p' Q {z) 



oo 
oo 

2 „/ (J\2 



^equimolar = 2lTl jdz(z- Z e ) p Q {z) . (44) 

— oo 

To evaluate these expressions, the density profiles po(z) and p\{z) still need to be 
determined from the expanded Euler-Lagrange equation: 

/'Oo) = /Wx + 2m pl(z) , (45) 

f'(po) Pl (z) = fn + 2mp'[{z) + 4mp' (z) . (46) 

In order to solve these equations, it is useful to assume proximity to the critical point 
so that the free energy density may be approximated by the usual double-well form: 

nrn 
f(p) ~ /WxP + Pcoex = , A ^ 2 (P - PO/) 2 (P - P0,v) 2 , (47) 

where the bulk correlation length £ is related to the second derivative of f(p) evaluated 
at either bulk density. Solving the Euler-Lagrange equation in Eq. (l4"5]) then leads to 
the usual tanh-form for the planar density profile [4j: 

p {z) = l -{ Po/ + p , v ) - ~y tanh((0 - z e )/2£) . (48) 

One may verify that solving the Euler-Lagrange equation in Eq. (146 j) gives the follow- 
ing general solution for pi(z) [15] : 

Pl (z) = ^m(Ap) 2 Z + ap' (z). (49) 

As already discussed, the rigidity constant is independent of the integration constant 
a. Inserting these profiles into the expressions for a, k and k in Eqs. (l4ip and 
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one finds fl 

m(Ap) 2 
a= ^^' (50) 

^equimolar = --(tT 2 - 3) 771 (Ap) 2 f , 
fcequimolar = -(7T 2 - 6) 771 (Ap) 2 £ . 

y 

For the symmetric double- well form for f(p), the Tolman length is identically zero. To 
obtain an estimate for 5 it is therefore necessary to consider leading order corrections 



to the double- well form for f(p) in E q.fETl) 14J . |34| . This leads to the following 
(constant) value for the Tolman length [34|: 



5 = -0.286565 ^m/a . (51) 

The prefactor depends on the precise form for f(p) and the number quoted is specific 
to the Carnahan-Starling equation of state 55]. 

All these formulas are derived assuming proximity to the critical point, but it turns 
out that they also provide a good approximation in a wide temperature range when 
the value of £ is chosen judiciously. This is done by using the fact that in squared- 
gradient theory the surface tension a may be determined from f(p) directly without 



fl: 



the necessity to determine the density profile po(z) 

Po,e 
2 y/m / dp \J f{p) - Pcoex P + Pcoex • (52) 



Po,e 
a ~\' 

P0,v 



An effective value for £ may now be chosen such that the two expressions for the 
surface tension in Eqs.f l50l) and ( |52l) are equal. This gives for £: 

£ — > £cff = — • 53 

3cr 

with a given by Eq.f l52|) . 

The procedure to determine the solid lines in Figures H] and [5] is now as follows. For 
a certain interaction potential, such as the Lennard- Jones potential in Eq.f f23|) . the 
interaction parameters a and m are calculated. Next, as a function of temperature, 
the bulk thermodynamic variables po,i, Po,v, Pcoex and p CO ex are derived from solving 
the set of equations in Eq. ([20l . The surface tension is then calculated from Eg. (152 j) 
and £ from Eq. (l53|) . With all parameters known, the curvature coefficients are finally 
calculated from Eqs. fl50|) and ( 15T|) . 
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V. LONG-RANGED INTERACTIONS: DISPERSION FORCES 

The surface tension, Tolman length and rigidity constants have all been explicitly 
evaluated using a Lennard- Jones potential that is cut-off beyond a certain distance r c . 
In this section we address the consequences of using the full Lennard- Jones potential. 
It is easily verified that the phase diagram in Figure [1] remains essentially the same 
when the cut-off is changed from r c = 7.5 to r c = oo, but that the shift in surface 
tension and Tolman length is increasingly noticeable (see Figured]). An inspection 
of the explicit expressions for the rigidity constants in Eq s.([29l) and fl30l) teaches us 



that both k and k diverge when r c increases to infinity [42], |51[ . This divergence is an 
indication that the expansion of the free energy is no longer of the form in Eq. ([6]) or 
(I7j), and it has to be replaced by 

25a , , T N log( d/R) , . 

a s (R) = a - — + (2k s + k s ) °y ; + . . . (54) 

v c (R) = a - — + k s ^ +... (55) 

where the dots represent terms of 0(1/ R 2 ). The coefficients of the logarithmic terms 
may be extracted from the expressions for k and k in Eqs.f l29l) and ( |30l) . They depend 
on the tail of the interaction potential, but are otherwise quite universal: 

k s = ^ed e (Ap) 2 , (56) 

h = -^edS(Ap) 2 . (57) 

This expression for k s is equal to that obtained in a DFT analysis of the singular part 



of the wave vector dependent surface tension of the fluctuating interface 56] . These 



expressions can also be derived from virial expressions for the rigidity constants when 



a sharp- kink approximation [5l| is made for the density profile [57] . The form for 



2k s + k s obtained by combining Eqs. (l56j) and (157"|) was first derived by Hooper and 



Nordholm in ref. 50] . 



To demonstrate the divergence of the second order term in Eq. (|54j) . the surface 
tension of a spherical liquid droplet as a function of the radius is determined for 
three values of the reduced LJ cut-off radius r c = 2.5, 7.5 and r c = oo. The regular 
contributions to <J S (R) from a and 5 are subtracted, so that we may define 

(2k + k) (R) = (a s (R) - a) R 2 + 25a R. (58) 
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FIG. 6: The combination (2k + k)(R) (in units of k-gT) as defined by Eq. (|58p as a function 
of the reciprocal equimolar radius d/R. The symbols are the results of DFT calculations at 
reduced temperature T* = 1.0 and three values of the reduced LJ cut-off r c = 2.5, 7.5 and 
oo. Solid circles are the corresponding values for 2k + k calculated from Eqs. (|34p and (|37p . 
The dashed line is the curve 7r/(6T*)(A/9*) 2 log(R /R) with i? - 0.005 d. 

This quantity is denned such that when the expansion in Eq.(j6]) for short-ranged 
forces is inserted, it reduces to 2k + k in the limit that R — > oo. For long-ranged 
forces (r c = oo), insertion of Eq.( |54l) into Eq.( |58l) gives a logarithmic divergence in 
this limit. This is verified by the DFT calculations shown in Figure [6] as the various 
symbols. For r c = 2.5 and r c = 7.5, the results indeed tend to the values obtained 
from the direct evaluation of 2 k + k using Eqs.( 1M|) and ( 137|) (solid circles). For 
r c = oo (triangular symbols) a slight divergence can be made out. This divergence is 
consistent with the dashed line, which is the divergence as described by combining 
the coefficients in Eqs.(l56l) and ( 1571) . 

VI. DISCUSSION 



In the context of density functional theory, we have shown that the surface tension 
of a spherical liquid droplet as a function of its inverse radius is well-represented 
by a parabola with its second derivative related to the rigidity constants k and k. 
Compact formulas for the evaluation of k and k are derived in terms of the density 
profiles po(z) and p\(z), which are in line with previous formulas presented by us 
34j and by Barrett [45|. A number of conclusions can be made with regard to these 
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formulas: 



The rigidity constants k and k depend on the choice for the location of the 
dividing surface of the planar density profile po(z). This dependency reflects 
the fact that when the location of the radius R is chosen differently, the curve of 
cr s (R) versus 1/R changes somewhat and the second derivative (2k+k) naturally 
needs to be amended. 

The most natural choice for a one-component system, is to locate the dividing 
surface of the planar interface according to the equimolar surface. For this choice 
both k and k are the least sensitive to a change in the location of the dividing 
surface. Furthermore, the equimolar value for k corresponds to its maximum 
value and the equimolar value for k corresponds to its minimum value. 

The bending rigidity k depends on the density profile pi(z), which measures the 
extent by which molecules rearrange themselves when the interface is curved. 
The bending rigidity is, however, independent of the choice made for the location 
of the dividing surface of p\ (z) (value of a in Eq. (136]) 



Using a cut-off and shifted Lennard- Jones potential for the attractive part of the 
interaction potential, the Tolman length and rigidity constants have been calculated 
with the result that S is negative with a value of minus 0.1-0.2 d, k is also negative with 
a value around minus 0.5-1.0 k-gT, and k is positive with a value of a bit more than 
half the magnitude of k. It is not expected that these results depend sensitively on the 
type of density functional theory used and we have shown that even an approximation 
scheme based on squared-gradient theory is quantitatively accurate. 

Our DFT results are expected to give an accurate qualitative description of the 
rigidity constants determined in experiments or computer simulations. First results 



of computer simulations by the group of Binder |24| shown in Figure El seem to sup- 
port this expectation, but further computer simulations are necessary. The agreement 
should cease to exist close to the critical point, however. Since the DFT calculations 
are all mean-field in character, the critical exponents obtained for both rigidity con- 
stants are the mean-field values of 1/2, which indicates that both k and k are zero at 
T c . Although it has not been proved rigorously, one expects that in reality the rigidity 
constants are finite at the critical point k, kock B T c . The situation is somewhat more 
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subtle for the rigidity constant associated with the description of surface fluctuations. 
Then, the bending rigidity is again negative but it vanishes on approach to the critical 



point with the same exponent as the surface tension 56]. 

The inspection of the explicit expressions presented for the rigidity constants is the 
most convincing method to investigate the possible presence of logarithmic correc- 



tions 47H49|. to replace the rigidity constants. For short-ranged interactions between 
molecules, the rigidity constants are definitely finite, but for an interaction potential 
that falls of as 1/r 6 for large intermolecular distances (dispersion forces), the rigidity 
constants are infinite indicating that the 1/R 2 term in the expansion of the surface 
tension needs to be replaced by a logarithmic term proportional to log(-R)/i? 2 . The 
proportionality constants of the logarithmic corrections are found to be quite univer- 
sal since they probe the systems long-distance behaviour and are in agreement with 



previous analyses [42 



50|,l5l|,|57|. 
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Appendix A: Alternative DFT expressions 

It may be useful to re-express the curvature coefficients 5, k, and k such that any 
reference to the chemical potential is absent. For the Tolman length the expression 
for Hi in Eq.( |25l) may be used to rewrite Eq. OTI) as: 

1 °° 
So- = _ Jd Zl Jdf l2 U M (r) r 2 (l - s 2 ) p^iK^) . (Al) 

— oo 

This expression is quite useful since it can be used to verify that the density profile 
Pi{z) determined numerically by solving the differential equation in Eq.( l25|) . leads to 
the same value for the Tolman length when evaluated using Eq. fl27j) . 

In order to transform the rigidity constants in a similar manner, we first need 
to expand the Euler-Lagrange equation in Eq. ffT9|) to second order in 1/R. For the 
spherical interface, one finds: 

/is, 2 = fL(Po) PsM) + \CiPo) Pi(^i) 2 + Jdri2 U att (r) [ Ps M) (A2) 
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+ y(l - s 2 )p[(z 2 ) - j(l-s 2 )z 2 p' (z 2 ) + j(l - s 2 y P >>{z 2 )] . 
The analogous expansion for the cylindrical interface gives: 

Pc,2 = fh s (po) Pc,2( z i) + -^fhs(po)pi(zi) 2 + Jdf 12 U att (r)[p C;2 (z 2 ) (A3) 

+ -(1 - S 2 ) P [(Z 2 ) - -(I - S 2 )Z 2P ' (Z 2 ) + —(1 - S 2 ) 2 Po(*2)] . 

Inserting these expressions for p s ^ 2 and p Cy2 into Eqs. (l29j) and (|30|) . one finds: 

oo 

k = - Jd Zl Jdf 12 U att (r) r 2 (l - s 2 ) p[,(zi)p^(^) (A4) 

— oo 
oo 

-^ Jd Zl Jdf 12 U att (r) r 2 (l - s 2 ) p' 1 (^i)p , 1 (z 2 ) 



-oo 
oo 



( /cfei /dri2 U att (r) r 2 (l - s 2 ) 2 z 2 p' (zi) p' (z 2 ) 



-oo 
oo 



+ ^ jdz 1 Jdf 12 £/ att (r)r 4 (l - s 2 )(l + 3s 2 ) p'^) p' Q {z 2 ) 



— oo 
oo 



fc = I Id: ., /^r 12 C/ a tt(r) r 2 (l - S 2 ) p(,( :,)[!//.,( V) -//,..,( : 2 )] (AO) 



-oo 
oo 



+ - /cfe x /dri2 ^att(^) ^ 2 (1 ~ s 2 ) «i Po(2i)po(z 2 ) 



-oo 
oo 



— oo 

These expressions have the advantage that no reference is made to the external field 
used to change the curvature. It might therefore be expected that these expressions 
are independent of the way the interfacial curvature is varied. An important disad- 
vantage, however, is that these expressions can only be evaluated when the second 
order corrections to the density profiles, p s , 2 (z) and p c ,2(z), are determined as well. 
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